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We study partially occupied lattice systems of classical magnetic dipoles which point along ran- 
domly oriented axes. Only dipolar interactions are taken into account. The aim of the model is to 
mimic collective effects in disordered assemblies of magnetic nanoparticles. From tempered Monte 
Carlo simulations, we obtain the following equilibrium results. The zero temperature entropy ap- 
proximately vanishes. Below a temperature Tc, given by fcsTc = (0.95 ± 0.1)a;eti, where Ed is a 
nearest neighbor dipole-dipole interaction energy and x is the site occupancy rate, we find a spin 
glass phase. In it, (1) the mean value (| q |), where q is the spin overlap, decreases algebraically with 
system size N as N increases, and (2) <5 | g |~ 0.5(| q \)^/tJx, independently of A'', where 5 | g | is 
the root mean square deviation of | g |. 
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I. INTRODUCTION 

Magnetic dipole-dipole interactions play a fundamen- 
tal role in magnetic phenomena. Their long range nature 
gives rise to magnetic domains in ferromagnets. On the 
other hand, because of their weak strength, they usually 
play an insignificant role in determining the Curie tem- 
peratures and in critical phenomena of atomic crystals. 
In assemblies of magnetic nanoparticles (NPs), things 
can be very different. Because ferromagnetic NPs of up 
to thousands of Bohr magnetons have a single magnetic 
domain,^ dipole-dipole interactions among such NPs can 
be very large, and can thus dominate their collective 
behavior.^ Often, both the orientation of the crystal- 
lites of which NPs are made as well as their positions 
are disordered in the assembly, and therefore behave 
much as a system of interacting magnetic dipoles with 
randomly oriented magnetic easy axes. These systems 
of random axes dipoles (RADs) are clearly frustrated, 
since two different dipoles give rise to magnetic fields at 
any given point which are not in general collinear. The 
sort of time dependent behavior that is expected of spin 
glasses has been observed in experimenta^'^'^ as well as 
in simulations^iii2ii2iii£ of disordered assemblies of mag- 
netic NPs. Because these systems evolve in time very 
slowly (exhibiting aging'^'- and other memory^ ^ effects), 
evidence for a thermodynamic spin glass phase in them 
is more difficult to come by. One of us has recently given 
numerical evidence for the existence of an equilibrium 
spin glass phase for a fully occupied lattice of dipoles 
with randomly oriented axes.^S, On the other hand, nu- 
merical evidence against such a phase has been given for 
site diluted lattice systems of magnetic dipoles with par- 
allel axesJ^^ 

One might expect a site diluted systems of RADs to be- 
have as fully occupied ones, and therefore to have equilib- 
rium spin glass phases at low temperatures. For support 



of this expectation, consider rescaling distances in a di- 
lute system of RADs. Because dipole dipole interactions 
decay with distance r as , letting r ^ br merely re- 
defines dipole-dipole interaction energies as Ed b~^ed- 
This would imply all physical quantities for systems of 
RAD's in three dimensions (3D) with different values of 
site concentration x collapse onto the same curve when 
plotted versus T/x, where T is the temperature. This 
argument holds for randomly located RAD's on a con- 
tinuous space in 3D. It is therefore expected to hold ap- 
proximately for a site diluted lattice if x ^ 1 , but not nec- 
essarily for lattices with higher concentrations of RADs. 
This is why dipoles with parallel axes can have an antifer- 
romagnetic phase at low temperature on a fully occupied 
simple cubic lattice)^ but appear to have no condensed 
phase of any sort on a very dilute lattice.— Similarly, the 
existence of a spin glass phase for RADs for a; < 1 does 
not follow from its existence for x = 1. 

Our main aim here is to find, by means of the parallel 
tempered Monte Carlo (TMC) algorithm^ii^ whether 
an equilibrium spin glass phase exists in a site diluted 
system of RADs in 3D. We also aim to establish, if the 
spin glass phase does exist, whether in the condensed 
phase (1) there is a single extended state^ii as in a ferro- 
magnet or in the droplet model of spin glasseSfi^ or (2) 
there are multiple extended states, as in the XY model in 
2D^- or in the replica symmetry breaking (RSB) theory 
of spin glasses 3 

The paper is planned as follows. In Sect. |TT]we specify 
the RAD model and describe how we apply the parallel 
tempered Monte Carlo (TMC) algorithmJ^ii^ We do this 
for X = 0.35 and x = 0.5. Results are given in Sect. IIIII 
The entropy S follows from our data for the specific heat 
by numerical integration. It approximately vanishes at 
zero temperature^^ More precisely, S < 0.01/cb at T = 0. 
We provide numerical evidence for the existence of an 
equilibrium spin glass phase below a transition temper- 
ature Tc- The evidence comes from the behavior of the 
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distribution of the spin overlap parameter q^^^ In anal- 
ogy to the behavior of the XY model in 2D, ( q^) seems 
to vanish in the macroscopic limit in the spin-glass phase, 
but only as a power of system size. For both x = 0.35 
and X = 0.5, fcsTc = (0.95 ± 0.1)xSd, where Sd is a near- 
est neighbor (NN) dipole-dipole interaction energy which 
is defined in Sect. |TT1 (Within errors, the value of for 
a; = l^'^ we have previously obtained also satisfies this 
expression.) For T < 0.9Tc, (5 | g |, that is, the root mean 
square deviation of | g |, fulfills S \ q \~ 0.5(| q \).JtJT^, 
independently of system size. This result is compared 
with its counterpart for the XY model in 2D. Results are 
discussed in Sect. Hvl 



II. MODEL AND METHOD 

We treat systems of magnetic dipoles on simple cu- 
bic lattices. We place a dipole on each lattice site with 
probability x, and leave the site unoccupied with prob- 
ability \ — X. Each dipole points along a randomly ori- 
ented anisotropy axes. With this model we aim to mimic 
assemblies of NPs in which uniaxial anisotropy energies 
are an order of magnitude larger than the largest dipole- 
dipole energy. Then, anisotropy energy barriers are not 
so large as to freeze spin reorientations near the spin glass 
temperature, but are sufficiently large to restrict spins to 
point "up" or "down" approximately along the easy mag- 
netization axes. We term this the random-axes-dipolar 
(RAD) model. The Hamiltonian is given by, 

ij af3 

where the first sum is over all occupied sites i and j of a 
simple cubic lattice, Sf is the a component of a classical 
3-component spin on site i, 

Tt/ = ed{a/n, f{5^p - -K^lrl;). (2) 

Tij is the distance between i and j, ea is an energy, and 
a a nearest neighbor distance. Each spin points along 
a randomly chosen direction. These equations can be 
cast into a form that is manifestly Ising like by letting 
Uj be (1) a null vector if the j site is unoccupied and 
(2) a 3— component unit vector chosen randomly from a 
spherically uniform distribution if the j site is occupied, 
and defining a pseudospin cjj = ±1 for each site, such 
that Sj = UjCTj. We can then write, 

'H = -^^Jij'^icrj, (3) 

where = -J2a,0T"fufu'^. Thus, the RAD model 
is an Ising model whose bonds are determined by 
the dipole-dipole terms T^"^ and the set of 3-component 
randomly oriented unit vectors {uj}. 



TABLE I: Simulation parameters, x is the probability that 
any given site is occupied by a magnetic dipole; A'^ is the mean 
number of magnetic dipoles in the system; AT is the temper- 
ature step in the TMC runs; Tm is the highest temperature 
of all systems; Ns is the number of disordered system pairs; 
MCS is the number of tempered Monte Carlo sweeps. 
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0.35 


0.35 


0.35 


0.35 


0.5 


0.5 


0.5 


0.5 


N 


22.4 


75.6 


179.2 


604.8 


32 


108 


256 


864 


AT 


0.05 


0.05 


0.05 


0.05 


0.05 


0.05 


0.05 


0.05 




1.0 


1.0 


1.0 


1.0 


1.3 


1.3 


1.3 


1.3 


Ns 


10000 


3000 


1250 


900 


5000 


3000 


1500 


150,600 


MCS 


50000 


10^ 


5 X lO'' 


lO'' 


50000 


50000 


10^ 


10^10^ 



From here on, unless we state otherwise, we let /c^ = 1, 
where fc^ is Boltzmann's constant, and give all temper- 
atures in terms of Ed- 

We use periodic boundary conditions. Details and jus- 
tification are given in Refs. [Uliiliillil . 

In order to arrive at equilibrium results, we make use of 
the parallel tempered Monte Carlo (TMC) algorithmfi^ 
This enables one to circumvent large energy barriers that 
can trap a system's state. We apply the TMC algorithm 
as follows. We run (in parallel) several identical systems 
at different temperatures: a system at temperature Tq, a 
second one at Tq + AT, and so on, at equally spaced tem- 
peratures, up to Tjn We choose Tm to be at least twice as 
large as what we expect to be the transition temperature 
between the paramagnetic and spin glass phases. We let 
each system evolve under the Metropolis MC algorithm 
for 10 MC sweeps before pairs of systems are given a 
chance to exchange their states. More specifically, pairs 
of systems at temperatures (Tq and Tq + AT), (To + 2 AT 
and To + SAT), and so on, are given a chance to ex- 
change states every 10 MC sweeps, and pairs of systems 
at temperatures {Tq + AT and To 2 AT), (To -h 3 AT and 
To -t- 4AT), and so on, are given a chance to exchange 
states at times in between. Under TMC rules, these 
chances are as follows. A system at temperature T that is 
in state 2, and another one at T-|- AT in state 1, exchange 
states with probability P = 1 if Ei < E2 and AT > 0, 
but P = exp{~Af3SE), where A/? = 1/T - 1/(T H- AT) 
and SE = El- E2, if Ei > E2 and AT > O.^^ 

A sufficiently small value of AT must be chosen in 
order to keep exp(— A/3(5i?) from becoming too small. 
This will often be fulfilled if A^AT; < 1, where AE 
is the mean energy difference between two systems at 
temperatures T and T -I- AT. The required condition, 
AT < T/ V A^c, where iV is the number of dipoles in the 
system, follows for AT if we replace AE by cN AT (we 
thus define c as the specific heat per spin). From plots 
of the specific heat vs T, such as the one shown in Fig. 
[l]for systems of 179 and of 605 magnetic dipoles on lat- 
tices with 0.35 of their sites occupied, one can get upper 
bounds for AT. (How the data points shown in Fig. [1] 
were obtained is explained in Sect. Illlh Values of AT as 
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FIG. 1: (a) (Color online) Plots of c and of 5" vs T, for x = 
0.35. Data points (black) o (for A*' = 605) and (red) x (for 
= 179) are for the specific heat, while (black) □ (for = 
605) and (red) + (for A'^ = 179) are for the entropy. Errors, 
5S, for S increase as T decreases, from SS ~ 0.002 for T = 4 
to 5S ^ 0.004 for T < 0.1. 



well as of other parameters for all TMC runs are given 
in table I. 

In order to probe for spin glass behavior, we define, as 
is usual, the spin overlap parameter .^^'^"^ First, let 



(4) 



where crj^-* and a^p are the pseudospins [defined above 
Eq. ([3])] on site j of identical twin replicas 1 and 2 of the 
system. Clearly, 



N- 



(5) 



is a measure of the spin configuration overlap between 
replicas 1 and 2. Thus, | (7 |= 1 if either crj^'' = crj^'' 
for all j or crj^'' ~ ^'^f^ We also define the 

moments of qk — {\ q \^), for fc = 1, 2 and 4, where 
(...) stands for an average which we next specify. 

Suppose our TMC runs last for a time t, and assume 
that, within the temperature range we are interested in, 
equilibration takes place in a time less than t/2. During 
the TMC runs, we write, at equally spaced time intervals 
within the (t/2,t) range, the spin configuration for each 
temperature. We later draw from these written configu- 
rations the numbers for that go into Eq. ([4]). This 
procedure is repeated for the same system (that is, for 
a system with the same anisotropy axes), but starting 
from different initial conditions. From this second set of 
configurations, we draw the numbers for a - that go into 
Eq. ^ . A value of q for each time t and temperature T 
is thus obtained from Eq. ([5]) . We first average all quan- 
tities of interest over and finally repeat these pairs of 
TMC runs a number Ng of times for different realizations 
of axes orientations, from which the average values we are 



reporting are obtained. Values we have used in our TMC 
runs for Ng and for t (labeled MCS therein) are given in 
table I. 

Finally, in order to check that equilibration actually 
takes place as assumed, we also calculate the spin overlap 
q, not between identical twin replicas, but between spin 
configurations at two different times tg and ti of the same 
TMC run. We do this for several values of tg and ti in the 
neighborhoods olt/2 and t respectively, and average over 
different random axes realizations, in order to obtain qi. 
For t sufficiently large, qi = qi- We find this is fulfilled 
for all parameter values given in table I, except for (1) 
X = 0.35, TV = 604.8 and T < 0.2, and (2) for x = 0.5, 
N = 864 and T < 0.35. 
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FIG. 2: (Color online) (a) Log-log plots of qi and of g2 vs A'^ 
for X — 0.35 and the shown values of T. Lines are guides to 
the eye. (b) Same as in (a) but for x = 0.5. Some data points 
for small T and large A'^ are missing, because our TMC runs 
did not reach equilibration for them. 



III. RESULTS 

Plots of the specific heat and of the entropy vs T for 
the RAD model are shown in Fig. [T]for x = 0.35. The 
data for the specific heat follows from numerical deriva- 
tives of the energy with respect to T. We obtain S 
from S{T) = ln2 + J^ciT')/T'dT' . Our data covers 
the temperature range (not all of it shown in Fig. [T]) 
0.05 < r < 4. For a numerical integration, we must ex- 
trapolate our data for c(T) beyond T = 4. To this end. 
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we use the leading term of an energy expansion in powers 
of 1/T, which gives c{T) A/T^. A fit of the value of 
A to our data for c{T) at T = 4 leads to a AS" ~ 0.027 
contribution to the entropy from the T > 4 range. We 
thus obtain S{T). It is exhibited in Fig. [TJ Errors for 
S{T) come in two approximately equal pieces: (1) an er- 
ror of roughly 0.002 from the 4 < T < oo range, and 
(2) an error of roughly 0.002, from errors in the data 
for the specific, which enter the integral c{T')/T'dT' . 
We finally obtain S = 0.015 ± 0.004 for T = 0.05, and 
extrapolations below T = 0.05 yield S < 0.01 at T = 0. 
Similarly, we find 5* < 0.01 at T = for x = 0.50. 
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FIG. 3: (Color online) (a) Plots of the probability distribution 
Pq vs Qr, for X — 0.35 and T = 0.3. •, o, and ♦, are for 
— 76, 179 and 605, respectively. Lines are guides to the eye. 
For clarity, only 30% of the data points are shown, (b) Same 
as in (a) but for T = 0.45. For comparison, a dashed line is 
shown for (I/tt) exp(— q^/tt), which ensues for a macroscopic 
paramagnetic phase. 

Our results for qi and q2 follow. Note we use an abso- 
lute value in the definition oiqi. Recall that Nq2 = 0{1) 
in the paramagnetic phase, and diverges at the transition 
temperature T^. Above the lower critical dimension dc, 
q2 = 0(1) in the droplet model of the spin glass phase. 
Our data shows that q2 decreases as N increases, for all 
nonzero temperatures. Log-log plots of q2 vs N are ex- 
hibited in Figs. [2^ and [2)d for various values of T. The 
behavior of q2 for T < 0.3 and T < 0.45 for x = 0.35 
and X = 0.5, respectively, is consistent with 92 ^ N^^, 
where y is some positive parameter that depends on T. 
This behavior is reminiscent of the XY model in 
For higher values of T, 52 vs A'^ clearly curves downwards, 
in accordance with a faster than algebraic in N decay, as 




FIG. 4: (Color online) (a) Plots of the probability Pq vs qr, 
for X = 0.50 and T = 0.45. •, o, and ♦, are for N = 108, 256 
and 864, respectively. Lines are guides to the eye. For clarity, 
only 30% of the data points are shown, (b) Same as in (a) 
but for T — 0.60. For comparison, a dashed line is shown for 
(I/tt) exp(— gj?/7r), which ensues for a macroscopic paramag- 
netic phase. 



one expects for the paramagnetic phase. Log-log plots 
of qf are also shown in Fig. [2^ and [2)3 for various val- 
ues of T. Note that (72/91 > 1 for all N and nonzero T, 
which implies a nonvanishing uncertainty in | g,. | , where 
Qr- = q/qi, for T > 0. 

We next report results for the probability distribution 
Pq{qr). Unless stated otherwise, all data given below 
are for a normalized Pq(qr), that is, J Pq{qr)dqr = 1. 
Recall that in the paramagnetic phase, because spin- 
spin correlation lengths are finite, the central limit the- 
orem implies qr is normally distributed for macroscopic 
systems. Plots of Pq vs qr are shown in Fig. 0^ (for 
T = 0.45) and Fig. (for T = 0.60), both for x = 0.50. 
The distribution of q appears to be size independent for 
T = 0.45. For T = 0.60, on the other hand, Pq{qr) 
drifts with system size. The drift seems consistent with 
Pqilr) (I/ti") exp(— (7,^/7r) as A*" ^ 00, which would 
be in accordance with a paramagnetic phase. Simi- 
lar remarks apply to the plots exhibited in Fig. [3^ 
(for T = 0.30) and Fig. Eb (for T = 0.45), both for 
X = 0.35. Clearly, 0.30 < < 0.45 for x ^ 0.35, and 
0.45 <Tc < 0.60 for x = 0.50. For T < T^, plots of Pq 
vs qr (not shown) become sharply peaked only as T — > 0. 
Now, if it turns out that dc = 3 for systems of RADs, then 
Pq{qr) ought to exhibit critical behavior for all T < Tc, 
more specifically, Pq{qr) ought to be size independent. 
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This sort of behavior is best summarized by the mean 
square deviation of | Qr \, = (qf) — (| |)^- 




FIG. 5: (Color online) (a) Plots of vs T for a; = 0.35. (red) 
• , (blue) ♦, and (black) □ are for N = 604.8, 179.2, and 75.6, 
respectively. The continuous (green) line is for A'^ = 22.4. 
The dashed line is for macroscopic paramagnet. (b) Same as 
in (a) but for x — 0.50. • (red), ♦ (blue), and □ (black) are for 
N = 864, 256, and 108, respectively. The continuous (green) 
line is for N = 32, and the dashed line is for 7r/2 — 1. (c) Same 
as in (a), but for the mean square relative deviation A\f of 
I M I, where M is the magnetic moment of XY systems of 
Lx L spins in 2D vs T/x, where x = 1, for L = 4 (green line), 
L = 8 (blue ♦), L = 16 (red •), L = 64 (black A), L = 256 
(black O), and L=1024 (black H). 

Plots of vs T/x are shown in Fig. (for x — 0.35) 
and in Fig. for x — 0.50 for various values of N. 
Clearly, curves in Figs. [5^ and [5Jd differ only slightly, 
as expected from the argument given in Sect. I, that 
variations of T and x have an effect only through T/x if 
x ^ 1. The data points in both figures suggest A^ 
7t/2 — 1 as — > oo, for T/x > 1. This is expected 
for a macroscopic paramagnet, since Pq{qr) is a normal 
distribution then. For T/x < 0.8, 

A, ~ 0.5^/t/^ (6) 



provides the best fit to the data points shown in Figs. [5^ 
and [5)3. This is in contrast with the behavior, A^ ^ 
as ^ 00, that one expects for the ordered phase of 
ferromagnets such as the Ising model in two or higher 
dimensions. For a more relevant comparison, we show in 
Fig. [SJ: plots of the mean square relative deviation A^ 
of I M I, where M is the magnetic moment, of the XY 
model in 2D. Notice there is no counterpart in Figs. [5^ 
and[5)D for the data points shown for rather large systems 
in Fig. [St. We note that A^ seems to go to 4 x 10~^T^-^ 
as TV ^ 00, for r < Tc in the XY model in 2D. 

Now, it is hard to see in Figs. [5^ and[S)3 where curves 
merge or cross, near T/x = 0.9. In order to enhance the 
differences between A^ curves for different values of L, we 
plot A^/A^(L = 4) vs T/x. These plots are shown in Fig. 
[S]for a; = 0.35 and x = 0.50. The weak dependence on 
X is remarkable. Notice that whereas curves for L ~ 12 
and L = 8 cross near T/x = 0.95 for x = 0.35 and for 
X = 0.50, the same curves cross the A^/A^(L = 4) = 1 
horizontal line only, if at all, for much smaller values of 
T/x. It thus appears that curves for A^ increasingly 
larger system sizes merge or cross at increasingly larger 
values of T/x. 




FIG. 6: (Color online) Plots of A^/A^(L = 4) vs T, where 
A^(L = 4) is A^ for i = 4. All (♦) are for x = 0.35 and 
L = 12, and aU (■) are for x = 0.35 and L = 8. AU (0) are 
for X = 0.50 and L — 12, and all (□) are for x — 0.50 and 
L = 8. Lines are guides to the eye. 

We show plots of ^^(0) vs A^ for a; = 0.35 and various 
temperatures in Fig. [71 In the paramagnetic phase, we 
expect Pq{Q) — > l/vr as A^ — > 00, as follows from a nor- 
mal distribution of Qr about qr = 0. Now, for T/x > 1.0, 
we see that Pq{0) < I/tt, but Pq{0) increases as A^ in- 
creases. On the other hand, for T/x < 0.857, Pq{0) 
remains, within errors, constant. This is again in con- 
trast with the behavior, Pq(0) — > as A^ — > 00, that 
one expects for the ordered phase of ferromagnets such 
as the Ising model in two or higher dimensions, but is in 
accordance with dc — 3. 

Finally, taking into account all our observations above, 
we arrive at our best estimate for Tc'. Tc/x = 0.95 ± 0.1. 
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IV. DISCUSSION 

By tempered Monte Carlo calculations, we have stud- 
ied systems of RADs in simple cubic lattices in which 
each site is occupied by a magnetic dipole with prob- 
ability X. Systems sizes, Monte Carlo run lengths and 
other details about the calculations can be found in Ta- 
ble I. The entropy S* as a function of temperature, which 
follows from our data for the specific heat, is shown in 
Fig. [1] for X — 0.35. S approximately vanishes at zero 
temperature. More precisely, S < O.Olfcs at T = for 
both X — 0.35 and x = 0.50j^ Note also that a vanishing 
zero temperature entropy implies that trapping above the 
lowest energy states as T — > O occurs only very rarely in 
our TMC runs. 

We have obtained the spin overlap parameter, defined 
in Eq. More specifically, we have obtained the mean 
value of I q I, which we term qi, and the mean square 
value of q, that is 92- The plots in Fig. [2^ (for x ~ 
0.35) and in Fig. [2)3 (for x = 0.50) suggest both qi 
and (72 vanish for all nonzero T as N 00. As can be 
gathered from those plots, both qf ~ N~y and 92 ^ N~y, 
where y is some positive parameter that depends on T, 
for ksT/x < 0.9ed- This is reminiscent of the XY model 
in 2D and suggests the lower critical dimension dc of the 
RAD model is at or near 3D. 

We have also studied the probability distribution 
Pq{qr), where g,- = q/qi- Within errors, Pq{qr) is in- 
dependent of system size at T = Tc, as is illustrated in 
Figs. [3^ andSK for ksT/x ~ O.Oe^. We have also found 
(not shown) Pq{qr) to be independent of system size at 
lower temperatures. The results shown in Figs. [S^ and 
and in Eq. ([6]) , for the mean square deviation of | | , 
A^, are consistent with a Pq{qr) that is independent of 
system size for all ksT/x < O.Qe^- Again, this is as in 
the XY model in 2D (see Fig. (5};) and suggests dc — 3 
for the RAD model. 



FIG. 7: Plots oiPq{q = 0) vs A^, for x = 0.35, for T/x = 1.286 
(A), T/x = 1.143 ( H), T/x = 1.0 ( ffl), T/x = 0.857 ( •), 
T/x = 0.714 ( ■), T/x = 0.57 ( ♦), and T/x = 0.429, (▼). 



There is a fairly large quantitative difference between 
the otherwise analogous behavior of the condensed phases 
of the XY model in 2D and the RAD model in 3D. 
From Eq. ([5]) for the RAD system, and from ~ 
4 X lO^'^r^'^, which we drew from the plots shown in Fig. 
[5t for the XY model, Ag/AA/ ~ 10^"°■^ where t = T/Tc, 
follows. Not surprisingly then, Pm/M at Af ~ for the 
XY model^5 is orders of magnitude smaller than Pq{0) for 
the RAD model, when they are both at some T/Tc well 
below a value of 1 . Replica symmetry breaking would ac- 
count this large quantitative difference between the XY 
and the RAD models at or near their lower critical di- 
mensionality. 
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